function re = au1_vf(p)
% AU1_VF Simulate Au1 HH Neurons

% initialize the results structures

% ib allocate
re.ib = {};
re.ib.so = {};
re.ib.ad = {};
re.ib.ax = {};

re.ib.so.series = zeros(p.ib.n,1+p.ib.so.n_current_param,p.n_steps+1);
re.ib.ad.series = zeros(p.ib.n,1+p.ib.ad.n_current_param,p.n_steps+1);
re.ib.ax.series = zeros(p.ib.n,1+p.ib.ax.n_current_param,p.n_steps+1);
re.ib.ad.ib_synapses = zeros(p.ib.n,p.ib.n,p.n_steps+1);
re.ib.ad.lts_synapses = zeros(p.ib.n,p.lts.n,p.n_steps+1);

% ib init
for i_ib = 1:p.ib.n  
    re.ib.so.series(i_ib,:,1) = initIB_SO(p.ib.so.cell{i_ib}.V_init);
    re.ib.ad.series(i_ib,:,1) = initIB_AD(p.ib.ad.cell{i_ib}.V_init);
    re.ib.ax.series(i_ib,:,1) = initIB_AX(p.ib.ax.cell{i_ib}.V_init);
end

% lts allocation
re.lts = {};

re.lts.series = zeros(p.lts.n,1+p.lts.n_current_param,p.n_steps);
re.lts.ib_synapses = zeros(p.lts.n,p.ib.n,p.n_steps);
re.lts.lts_synapses = zeros(p.lts.n,p.lts.n,p.n_steps);

% lts init
for i_lts = 1:p.lts.n
    re.lts.series(i_lts,:,1) = initLTS(p.lts.cell{i_lts}.V_init);
end

% loop through time
for t_step = 2:p.n_steps+1
    V_lts_all = re.lts.series(:,1,t_step-1);
    V_ib_ax_all = re.ib.ax.series(:,1,t_step-1);
    for i_lts = 1:p.lts.n
        lts_next = rk4_solver(@LTS,p.dt,p.time(t_step-1),[re.lts.series(i_lts,:,t_step-1),re.lts.lts_synapses(i_lts,:,t_step-1),re.lts.ib_synapses(i_lts,:,t_step-1)],p,i_lts,V_lts_all,V_ib_ax_all);
        re.lts.series(i_lts,:,t_step) = lts_next(1:1+p.lts.n_current_param);
        re.lts.lts_synapses(i_lts,:,t_step) = lts_next(1+p.lts.n_current_param+1:1+p.lts.n_current_param+p.lts.n);
        re.lts.ib_synapses(i_lts,:,t_step) = lts_next(1+p.lts.n_current_param+p.lts.n+1:end);
    end
    for i_ib = 1:p.ib.n
        V_ib_so = re.ib.so.series(i_ib,1,t_step-1);
        V_ib_ad = re.ib.ad.series(i_ib,1,t_step-1);
        V_ib_ax = re.ib.ax.series(i_ib,1,t_step-1);
            
        ib_ad_next = rk4_solver(@IB_AD,p.dt,p.time(t_step-1),[re.ib.ad.series(i_ib,:,t_step-1),re.ib.ad.lts_synapses(i_ib,:,t_step-1),re.ib.ad.ib_synapses(i_ib,:,t_step-1)],p,i_ib,V_ib_so,V_lts_all);
        re.ib.ad.series(i_ib,:,t_step) = ib_ad_next(1:1+p.ib.ad.n_current_param);
        re.ib.ad.lts_synapses(i_ib,:,t_step) = ib_ad_next(1+p.ib.ad.n_current_param+1:1+p.ib.ad.n_current_param+p.lts.n);
        re.ib.ad.ib_synapses(i_ib,:,t_step)= ib_ad_next(1+p.ib.ad.n_current_param+p.lts.n+1:end);
        
        re.ib.so.series(i_ib,:,t_step) = rk4_solver(@IB_SO,p.dt,p.time(t_step-1),re.ib.so.series(i_ib,:,t_step-1),p,i_ib,V_ib_ad,V_ib_ax);
        
        re.ib.ax.series(i_ib,:,t_step) = rk4_solver(@IB_AX,p.dt,p.time(t_step-1),re.ib.ax.series(i_ib,:,t_step-1),p,i_ib,V_ib_so);
    end
end

end

% UPDATE FUNCTIONS

% LTS Cells

function dy = LTS(t, y, p, i_lts, V_lts, V_ib_ax)
%LTS defines the dynamics of a superficial LTS interneuron

n_y_inputs = length(y);
if n_y_inputs ~= p.lts.n + p.ib.n + p.lts.n_current_param + 1;
    error('LTS: number of inputs does not match expected number of inputs');
end

V = y(1);
h_L = y(2);
m_L = y(3);
m_ARL = y(4);
s_lts = y(1+p.lts.n_current_param+1:1+p.lts.n_current_param+p.lts.n);
s_ib = y(1+p.lts.n_current_param+p.lts.n+1:end);

dy = zeros(1,n_y_inputs); 

% intrinsic currents
I_Leak = p.lts.g_Leak * (65 + V);
I_NaF = p.lts.g_NaFL * (m_0i(V)^3) * h_L * (V - 50);
I_KDR = p.lts.g_KDRL * (m_L^4) * (100 + V);
I_AR = p.lts.g_ARL * m_ARL * (35 + V);

% synaptic currents
Is_LTS = p.lts.cell{i_lts}.lts_gate .* s_lts .* (V + p.lts.Vk_L);
Is_IB = p.lts.cell{i_lts}.ib_gate .* s_ib .* (V + p.lts.Vk_a);

% dV
dy(1) = (-1/p.lts.C) * (p.lts.cell{i_lts}.J + I_Leak + I_NaF + I_KDR + I_AR + sum(Is_LTS) + sum(Is_IB));
% dh_L
dy(2) = (1/tau_hi(V))* (h_infi(V) - h_L);
% dm_L
dy(3) = (1/tau_mi(V)) * (m_infi(V) - m_L);
% dm_ARL
dy(4) = (1/tau_mAR(V)) * (m_infAR(V,75) - m_ARL);

% ds_lts
dy(1+p.lts.n_current_param+1:1+p.lts.n_current_param+p.lts.n) = ...
    synapse(s_lts, p.lts.tau_sLL_r, p.lts.tau_sLL_d, V_lts);
% ds_ib
dy(1+p.lts.n_current_param+p.lts.n+1:end) = ...
    synapse(s_ib, p.lts.tau_saL_r, p.lts.tau_saL_d, V_ib_ax);

end

% IB Cells - Apical Dendrite

function dy = IB_AD(t, y, p, i_ib, V_ib_so, V_lts)
%IB_AD defines the dynamics of the dendrite of a deep IB cell.

n_y_inputs = length(y);
if n_y_inputs ~= p.ib.n + p.lts.n + p.ib.ad.n_current_param + 1;
    error('IB_AD: number of inputs does not match expected number of inputs');
end

V = y(1);
h_d = y(2);
m_d = y(3);
m_ARd = y(4);
m_KMd = y(5);
m_CaHd = y(6);
s_lts = y(1+p.ib.ad.n_current_param+1:1+p.ib.ad.n_current_param+p.ib.n);
s_ib = y(1+p.ib.ad.n_current_param+p.ib.n+1:end);

dy = zeros(1,n_y_inputs); 

% instrinsic currents
I_Leak = p.ib.ad.g_Leak * (70 + V);
I_NaF = p.ib.ad.g_NaFd * (m_0e(V)^3) * h_d * (V - 50);
I_KDR = p.ib.ad.g_KDRd * (m_d^4) * (95 + V);
I_AR = p.ib.ad.g_ARd * m_ARd * (25 + V);
I_KM = p.ib.ad.g_KMd * m_KMd * (95 + V);
I_CaH = p.ib.ad.g_CaHd * (m_CaHd^2)*(V - 125);
Is_sd = -1 * p.ib.ad.g_sd * (V_ib_so - V); % soma to dendrite

% synaptic currents
Is_LTS = p.ib.ad.cell{i_ib}.lts_gate .* s_lts * (V + p.ib.ad.Vk_L);
Is_IB = p.ib.ad.cell{i_ib}.ib_gate .* s_ib * (V + p.ib.ad.Vk_a);

% dV
dy(1) = (-1/p.ib.ad.C) * (p.ib.ad.cell{i_ib}.J + I_Leak + I_NaF + I_KDR + I_AR + I_KM + I_CaH + Is_sd + sum(Is_LTS) + sum(Is_IB));
% dh_d 
dy(2) = (1/tau_he(V)) * (h_infe(V) - h_d);
% dm_d
dy(3) = (1/tau_me(V)) * (m_infe(V) - m_d);
% dm_ARd
dy(4) = (1/tau_mAR(V)) * (m_infAR(V,75) - m_ARd);
% dm_KMd
dy(5) = alpha_KM(V) * (1-m_KMd) - beta_KM(V) * m_KMd;
% dm_CaHd
dy(6) = (1/.33) * (alpha_CaH(V) * (1 - m_CaHd) - beta_CaH(V) * m_CaHd);
% ds_lts
dy(1+p.ib.ad.n_current_param+1:1+p.ib.ad.n_current_param+p.lts.n) = ...
    synapse(s_lts, p.ib.ad.tau_sLd_r, p.ib.ad.tau_sLd_d, V_lts);
% ds_ib - todo - cmk 5/17/10 - not implemented for ib to ib connections
dy(1+p.ib.ad.n_current_param+p.ib.n+1:end) = zeros(p.ib.n,1);

end

% IB Cells - Axon

function dy = IB_AX(t, y, p, i_ib, V_ib_so)
%IB_A defines the dynamics of the axon of a deep IB cell.

n_y_inputs = length(y);
if n_y_inputs ~= p.ib.ax.n_current_param + 1;
    error('IB_AX: number of inputs does not match expected number of inputs');
end

V = y(1);
h_a = y(2);
m_a = y(3);
m_KMa = y(4);

dy = zeros(1,n_y_inputs); 

I_Leak = p.ib.ax.g_Leak * (70 + V);
I_NaF = p.ib.ax.g_NaFa * (m_0e(V)^3) * h_a * (V-50);
I_KDR = p.ib.ax.g_KDRa * (m_a^4) * (95 + V);
I_KM = p.ib.ax.g_KMa * m_KMa * (95 + V);
Is_sa = -1 * p.ib.ax.g_sa * (V_ib_so - V); % soma to axon

% note - cmk 5/18/10 - removed Ic_gj - gap junction voltage not happening in this model

% dV 
dy(1) = (-1/p.ib.ax.C) * (p.ib.ax.cell{i_ib}.J + I_Leak + I_NaF + I_KDR + I_KM + Is_sa);
% dh_a
dy(2) = (1/tau_he(V)) * (h_infe(V) - h_a);
% dm_a
dy(3) = (1/tau_me(V)) * (m_infe(V) - m_a);
% dm_KMa
dy(4) = alpha_KM(V) * (1-m_KMa) - beta_KM(V) * m_KMa;

end

% IB Cells = Soma

function dy = IB_SO(t, y, p, i_ib, V_ib_ad, V_ib_ax)
%IB_BD defines the dynamics of the soma of a deep IB cell.

n_y_inputs = length(y);
if n_y_inputs ~= p.ib.so.n_current_param + 1;
    error('IB_SO: number of inputs does not match expected number of inputs');
end

V = y(1);
h_s = y(2);
m_s = y(3);

dy = zeros(1,n_y_inputs); 

I_Leak = p.ib.so.g_Leak * (70 + V);
I_NaF = p.ib.so.g_NaFs * (m_0e(V)^3) * h_s * (V-50);
I_KDR = p.ib.so.g_KDRs * (m_s^4) * (95 + V);
Is_as = -1 * p.ib.so.g_as * (V_ib_ax - V); % axon to soma
Is_ds = -1 * p.ib.so.g_ds * (V_ib_ad - V); % dendrite to soma

% dV
dy(1) = (-1/p.ib.so.C) * (p.ib.so.cell{i_ib}.J + I_Leak + I_NaF + I_KDR + Is_as + Is_ds);
% dh_s
dy(2) = (1/tau_he(V)) * (h_infe(V) - h_s);
% dm_s
dy(3) = (1/tau_me(V)) * (m_infe(V) - m_s);

end

% constant calculators for different currents

% m currents

function m = m_0i(V)
%m_0e calculates the rate constant m0 for excitatory cell NaF & KDR
%currents
%   Detailed explanation goes here

m = 1/(1+exp((-V-38)/10));

end

function m = m_infi(V)
%m_infe calculates the rate constant m_inf for excitatory cell NaF & KDR
%currents
%   Detailed explanation goes here

m = 1/(1+exp((-27-V)/11.5));

end

function m = m_infe(V)
%m_infe calculates the rate constant m_inf for excitatory cell NaF & KDR
%currents
%   Detailed explanation goes here

m = 1/(1+exp((-29.5-V)/10));

end

function m = m_infAR(V, V0)
%m_infAR calculates the rate constant m_infAR for h-currents
%   Detailed explanation goes here

m = 1/(1+exp((V+V0)/5.5));

end

function m = m_0e(V)
%m_0e calculates the rate constant m0 for excitatory cell NaF & KDR
%currents
%   Detailed explanation goes here

m = 1/(1+exp((-34.5-V)/10));

end

% h currents

function h = h_infi(V)
%h_infe calculates the rate constant h_inf for excitatory cell NaF & KDR
%currents
%   Detailed explanation goes here

h = 1/(1+exp((V+58.3)/6.7));

end

function h = h_infe(V)
%h_infe calculates the rate constant h_inf for excitatory cell NaF & KDR
%currents
%   Detailed explanation goes here

h = 1/(1+exp((V+59.4)/10.7));

end

function t = tau_mi(V)
%tau_he calculates the rate constant tau_h for excitatory cell NaF & KDR
%currents

t = 0.25+4.35*exp((-1*abs(V+10))/10);

end

function t = tau_me(V)
%tau_he calculates the rate constant tau_h for excitatory cell NaF & KDR
%currents

t = 0.25+4.35*exp((-1*abs(V+10))/10);

end

function m = tau_mAR(V)
%m_infAR calculates the rate constant m_infAR for h-currents
%   Detailed explanation goes here

m = 1/(exp(-14.6-.086*V)+exp(-1.87+.07*V));

end

function t = tau_hi(V)
%tau_he calculates the rate constant tau_h for excitatory cell NaF & KDR
%currents
%   Detailed explanation goes here

t = 0.225 + (1.125/(1+exp((V+37)/15)));

end

function t = tau_he(V)
%tau_he calculates the rate constant tau_h for excitatory cell NaF & KDR
%currents
%   Detailed explanation goes here

t = 0.15 + (1.15/(1+exp((V+33.5)/15)));

end

% synapse dynamics

function ds = synapse(s,tau_r,tau_d,V_pre)
%ds returns the dynamics of a synapse
%   input the s variable, time constant for rise and decay, and pre-
%   synaptic voltage.

ds = (-s/tau_d) + (((1-s)/tau_r) * (1+tanh(V_pre/10)));

end

% initial conditions

% lts init

function lts0 = initLTS(V_init)
%initLTS generates an array of numCells LTS cells at initial conditions

lts0 = [V_init h_infi(V_init) m_infi(V_init) m_infAR(V_init,75)];

end

% ib init

function ib_ad0 = initIB_AD(V_init)

m_KM_inf = alpha_KM(V_init) / (alpha_KM(V_init) + beta_KM(V_init));
m_CaH_inf = alpha_CaH(V_init) / (alpha_CaH(V_init) + beta_CaH(V_init));

ib_ad0 = [V_init h_infe(V_init) m_infe(V_init) m_infAR(V_init, 75) m_KM_inf m_CaH_inf];

end

function ib_so0 = initIB_SO(V_init)

ib_so0 = [V_init h_infe(V_init) m_infe(V_init)];

end

function ib_ax0 = initIB_AX(V_init)

m_KM_inf = alpha_KM(V_init)/(alpha_KM(V_init) + beta_KM(V_init));

ib_ax0 = [V_init h_infe(V_init) m_infe(V_init) m_KM_inf];

end

function a = alpha_KM(V)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

a = .02/(1.0+exp((-V-20)/5));

end

function a = alpha_CaH(V)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

a = 1.6/(1.0+exp(-.072*(V-5)));

end

function b = beta_KM(V)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

b = .01*exp((-V-43)/18);

end

function b = beta_CaH(V)
%UNTITLED Summary of this function goes here
%   Detailed explanation goes here

b = (.02*(V+8.9))/(exp((V+8.9)/5)-1);

end

% rk4_solver
% runs a vector of variables simultaneously through the runge-kutta 4
% method of solving differential equations

function y_next = rk4_solver(fn,h,t,y,varargin)
 
    k1 = (h * feval(fn, t, y, varargin{:}));
    k2 = (h * feval(fn, t + h/2, y + k1/2, varargin{:}));
    k3 = (h * feval(fn, t + h/2, y + k2/2, varargin{:}));
    k4 = (h * feval(fn, t + h, y + k3, varargin{:}));
    y_next = y + (k1 + 2*k2 + 2*k3 + k4)/6;

end